function [avgIm] = computeDirAvgImage(imgPath, fileTypeSpec, methodType)
%computeDirAvgImage
%   Performs pre-op work on images for algorithm. If input is a directory
%   of images the average of all images in folder is computed. If input is
%   a single image, the image is converted into 2D and double.
% 
%   Inputs
%   ---------
%   imgPath may either be the folder_name output from uigetdir or a
%   fileName output from uigetfile. It reads in images from a filename or
%   path name
%   
%   fileTypeSpec is the type of image files that the user wants to average.
%
%   methodType determines whether to use the single image or directory
%   protocols for image prep. methodType can take a string of 'file' or
%   'dir'.
%
%   Output
%   ----------
%   avgIm is the output image from the pre-op work. If input is an image file,
%   this will be a 2D image that is class double. If input is a directory
%   of images, avgIm will be the average of the images in the directory.

%Single Image Case
if strcmp(methodType, 'file')
    singleImage = imread(imgPath);
    if ndims(singleImage) == 3
        singleImage = rgb2gray(singleImage); %Remove 3rd dimension for compatibility
    end
    singleImage = im2double(singleImage); %Convert to double for algorithm compatibility
    avgIm = singleImage;
elseif strcmp(methodType, 'dir')
    imageFiles = dir([imgPath filesep fileTypeSpec]);
    num_images = length(imageFiles);
    if (num_images == 0)
        error('Please select a directory with registered images.')
    end
    imSize = size(imread([imgPath filesep imageFiles(1).name]));
    imSize = imSize([1 2]); % ignore 3rd dimension for rgb images.
    
    imCube = zeros(imSize(1),imSize(2),num_images);
    
    for i=1:num_images
        temp_im = imread([imgPath filesep imageFiles(i).name]);
        if ndims(temp_im) == 3
            temp_im = rgb2gray(temp_im);
        end
        temp_im = im2double(temp_im);
        
        if(size(temp_im) == imSize)
            imCube(:,:,i) = temp_im;
        else
            error('Please choose a directory with only registered images of the same dimensions.');
        end
    end
    % Compute overlapping mask so only average of valid pixel locations is taken
    count = 1;
    for i = 1:size(imCube, 3)
        imgbwCube(:,:,count) = im2bw(imCube(:,:,i), 0.09);
        count = count + 1;
    end
    for i = 2:num_images
        imgbwMask = immultiply(imgbwCube(:,:,1), imgbwCube(:,:,i));
    end
    
    % Multiply mask onto all images
    % Multiply by same class and need to have type double
    imgbwMask = im2double(imgbwMask);
    imCube = im2double(imCube);
  
    imgcropArray = zeros([size(imCube, 1), size(imCube, 2), num_images]);
    for k = 1:num_images
        imgcropArray(:,:,k) = imCube(:,:,k).*imgbwMask;
    end
    
    % Compute Average
    avgIm = mean(imgcropArray,3);
 
else
    error('Please check that a proper string for the methodType input was used')
end

end

